# The following script was adopted from Thiele, Jan C., Kurth, Winfried and Grimm, Volker (2014) 
# 'Facilitating Parameter Estimation and Sensitivity Analysis of Agent-Based Models: 
# A Cookbook Using NetLogo and 'R'' Journal of Artificial Societies and Social Simulation 17 (3) 11 
# <http://jasss.soc.surrey.ac.uk/17/3/11.html>.


library(RNetLogo)

# an R random seed (for beeing reproducible)
set.seed(-1589217668)

# NetLogo path; NetLogo model path; simulation function path

nl.path     <- "C:/Program Files (x86)/NetLogo 5.2.0"
model.path  <- "C:/.../ABSAM_R.nlogo"
simfun.path <- "C:/.../simulation_function.R"

# parameters to be calibrated
parameter.values <- list(
                          'shock'  = list(min = 0.01, max = 0.04, random.function="qunif"),
                          'minimum-wage' = list(min = 0.8 , max = 1.7, random.function="qunif"),
                          'efficiency'    = list(min = 0.10,  max = 0.30, random.function="qunif"),
                          'beta'          = list(min = 0.4, max = 0.6, random.function="qunif"),
                          'growth-rate'    = list(min = 0.002,  max = 0.03, random.function="qunif"),
                          'benefits'      = list(min = 0.4, max = 1.2, random.function="qunif")
                        )

# Min and max of calibration criteria extracted from empirical data
calibration.ranges <- data.frame(
                            'u_rate' = c(6.4 , 14.6),
                            'theta' = c(0.15 , 0.35),
                            'ltu_rate'   = c(19.7 , 47.2)
                          )

# number of parameter sets
sample.count <- 120

# how many repetitions for each parameter set should be run
no.repeated.sim <- 10

# should R report the progress
trace.progress = TRUE

# initialize NetLogo
nl.LMA2.0R.nlogo <- "nl.LMA2.0R.nlogo"
NLStart(nl.path, gui=TRUE, nl.obj=nl.LMA2.0R.nlogo)
NLLoadModel(model.path,nl.obj=nl.LMA2.0R.nlogo)

# load the code of the simulation function
source(file=simfun.path)

# get names of parameters
parameter.names <- names(parameter.values)

library(lhs)

# create a random sample of parameter sets
lhs.design <- randomLHS(sample.count, length(parameter.values))

# transform the standardized random values to the real parameter value range
lhs.design <- lapply(seq(1,length(parameter.values)), function(i) {
    match.fun(parameter.values[[i]]$random.function)(lhs.design[,i], parameter.values[[i]]$min, 
                                                     parameter.values[[i]]$max)
  })
names(lhs.design) <- parameter.names


already.processed <- 0

# simulate for all parameter sets and get results of all evalulation criteria as matrix
sim.results.lhs <- apply(as.data.frame(lhs.design), 1, 
                         simulate, no.repeated.sim=no.repeated.sim, 
                         nl.obj=nl.LMA2.0R.nlogo, trace.progress=trace.progress,
                         parameter.names=parameter.names, 
                         iter.length=sample.count,
                         function.name="LHS"
                        )

# get names of calibration criteria
calibration.crit.names <- names(calibration.ranges)

# check, if simulation results are in the calibration ranges for each parameter and simulation
eval.param <- data.frame(sim.results.lhs > t(calibration.ranges)[,1] & sim.results.lhs < t(calibration.ranges)[,2])
row.names(eval.param) <- calibration.crit.names

# how many criteria were met for each simulation?
count.criteria.met <- colSums(eval.param)
summary(count.criteria.met)

# is there a parameter set which fullfils all criteria?
all.criteria.met <- which(colSums(eval.param) == length(eval.param[,1]))
print(all.criteria.met)

# print the successfull parameter sets
if (length(all.criteria.met) > 0)
{
  as.data.frame(lhs.design)[all.criteria.met,]
}


# plot the fullfilled criteria

library(scatterplot3d)
plot.calibr.points <- function(critno) {
  matched <- which(t(eval.param)[,critno]==T)  
  mp$points3d(lhs.design[[1]][matched], lhs.design[[2]][matched], lhs.design[[3]][matched], pch=critno, col=1+critno, cex=4.5, lwd=3)
}

plot.calibr.points1 <- function(critno) {
  matched <- which(t(eval.param)[,critno]==T)  
  mp1$points3d(lhs.design[[4]][matched], lhs.design[[5]][matched], lhs.design[[6]][matched], pch=critno, col=1+critno, cex=4.5, lwd=3)
}

par(mfrow=c(2,1), mar=c(0,0,0,0), oma = c(0,0,0,0), mgp=c(5,1,3))
mp <- scatterplot3d(lhs.design[[1]], lhs.design[[2]], lhs.design[[3]], 
                    color="black", lwd=4, pch=15, col.axis="black",
                    xlab="shocks frequency", ylab=parameter.names[2], zlab=parameter.names[3],
                    y.margin.add=0.1, col.lab="black", lty.grid="dotted", col.grid="black",
                    box = TRUE, scale.y = 1, angle = 25, mar=c(3,3,0,3),
                    axis=TRUE, cex.axis=1.1, cex.lab=1.3, tick.marks=TRUE)

invisible(lapply(1:length(calibration.crit.names), plot.calibr.points))

mp1 <- scatterplot3d(lhs.design[[4]], lhs.design[[5]], lhs.design[[6]], 
                     color="black", lwd=4, pch=15, col.axis="black",
                     xlab=parameter.names[4], ylab=parameter.names[5], zlab=parameter.names[6],
                     y.margin.add=0.1, col.lab="black", lty.grid="dotted", col.grid="black",
                     box = TRUE, scale.y = 1, angle = 25, mar=c(3,3,0,3),
                     axis=TRUE, cex.axis=1.1, cex.lab=1.3, tick.marks=TRUE)

# plot criteria points
invisible(lapply(1:length(calibration.crit.names), plot.calibr.points1))
